Helium in the Extended Atmosphere of the Warm Superpuff TOI-1420b

Thirty years of exoplanet discovery have revealed planetary archetypes for which we have no direct counterparts in our solar system. Among these myriad surprises are the superpuffs, a class of low-mass planets with densities ρ ≲ 0.1 g cm⁻³ (Lee & Chiang 2016; Jontof-Hutter 2019). For planets with M_p ≲ 50M_⊕, such low densities imply core masses smaller than 5M_⊕, which are challenging to reconcile with the core accretion theory for planet formation unless the envelopes are accreted "dust-free" (Lee & Chiang 2016) or atop ice-rich cores (Stevenson 1984; Venturini et al. 2015) beyond the water ice line. Moreover, if superpuffs are indeed as gas-rich as their densities imply, then many should be prone to catastrophic mass loss via photoevaporation and/or boil-off, with atmospheric lifetimes shorter than or comparable to the system ages (i.e., ${M}_{p}/\dot{M}\lesssim$ 1 Gyr for catastrophic loss; Cubillos et al. 2017; Wang & Dai 2019; Chachan et al. 2020). A number of mechanisms have been proposed to help ease these tensions and bring superpuff envelope mass fractions into closer agreement with expectations from core accretion, including inflation by tides (Millholland 2019; Millholland et al. 2020) or heat from formation (Masuda 2014; Libby-Roberts et al. 2020), planetary rings (Akinsanmi et al. 2020; Piro & Vissapragada 2020; Ohno & Fortney 2022), and high-altitude aerosols (Kawashima et al. 2019; Wang & Dai 2019; Gao & Zhang 2020; Ohno & Tanaka 2021).

We show the superpuffs in context on an M_p − T_eq diagram in Figure 1. All transiting planets on the NASA Exoplanet Archive (Akeson et al. 2013) with densities ρ < 0.3 g cm⁻³ and masses measured to better than 30% are shown, and planets on grazing orbits have been removed. Most puffy planets hotter than 1000 K are typical gas giants, with large H₂/He envelopes that are readily explained by core accretion models (e.g., Lee 2019). Gas giants above this temperature threshold have long been known to have inflated radii driven by the large incident stellar flux (e.g., Bodenheimer et al. 2003; Thorngren & Fortney 2018). Planets usually considered "superpuffs" (i.e., those with inferred core masses too small to trigger classical runaway gas accretion) are found in the lower left of this figure. The superpuff and gas giant regimes appear to be separated by a deficit of planets, which may result from sensitivity biases as a function of the orbital period: most gas giants have masses measured with radial velocities, and most superpuffs with transit timing variations (Steffen 2016; Mills & Mazeh 2017).

Zoom In Zoom Out Reset image size

The population is shown in more detail in the right panel of Figure 1. We focus on the planets in the upper right corner of the plot, WASP-107b (Anderson et al. 2017) and TOI-1420b (Yoshida et al. 2023). Both have inferred core masses <5M_⊕ (Piaulet et al. 2021; Yoshida et al. 2023), too small to trigger runaway gas accretion in the traditional core accretion paradigm for gas giants (which requires core masses ≳10M_⊕; e.g., Bodenheimer & Pollack 1986). However, these planets stand out as warmer and higher-mass than the other superpuffs. They also differ from the rest of the superpuff population in system architecture: WASP-107b is accompanied only by a distant nonresonant companion (Piaulet et al. 2021), and TOI-1420b also lacks nearby companions (Yoshida et al. 2023). All other superpuffs are found in multiplanet systems near mean-motion resonances, which facilitates their characterization via transit timing variations (Jontof-Hutter 2019), although a few have now also been detected with radial velocities (Petigura et al. 2018; Santerne et al. 2019).

Although their short orbital periods (P < 10 days) make them more prone to XUV-driven photoevaporation, TOI-1420b and WASP-107b are not expected to be susceptible to catastrophic atmospheric loss, as their gravitational potentials are relatively large despite their low densities (Wang & Dai 2019; Chachan et al. 2020; Belkovski et al. 2022). Indeed, He 1083 nm measurements of WASP-107b suggest an outflow rate of about one Earth mass per Gyr ($M/\dot{M}\approx 30\,\mathrm{Gyr};$ Spake et al. 2018; Allart et al. 2019; Kirk et al. 2020; Spake et al. 2021; Wang & Dai 2021), in line with predictions from photoevaporation models (Kubyshkina et al. 2018; Caldiroli et al. 2022). This planet has not lost much envelope mass over its lifetime, so it remains plausible that WASP-107b formed with an unusually large envelope (Lee & Chiang 2016; Piaulet et al. 2021). This is a major difference from the lower-mass superpuffs, which would be expected to lose large envelopes on short timescales necessitating many of the aforementioned alternative mechanisms (rings, hazes, etc.).

In this work, we study atmospheric escape in TOI-1420b (Yoshida et al. 2023). This planet orbits a late G dwarf every 6.96 days, and with M_p = 25.1 ± 3.8M_⊕ and R_p = 11.9 ± 0.3R_⊕, TOI-1420b has an exceptionally low density of just 0.08 ± 0.02 g cm⁻³. Moreover, at J = 10.6 the system is well-positioned for narrowband photometric observations with the Wide-field InfraRed Camera (WIRC; Wilson et al. 2003) on the Hale 200-inch Telescope at Palomar Observatory. We recently used Palomar/WIRC to constrain helium absorption in three systems with similar J magnitudes: HAT-P-18 (J = 10.8), WASP-52 (J = 10.6), and HAT-P-26 (J = 10.1; Vissapragada et al. 2020b; Paragas et al. 2021; Vissapragada et al. 2022b). The former two measurements have recently been confirmed by JWST NIRISS/SOSS and Keck/NIRPSEC, demonstrating the utility of narrowband photometry for this science case (Fu et al. 2022; Kirk et al. 2022). In Section 2, we describe our narrowband photometric observations of TOI-1420b. We model our narrowband helium light curves in Section 3, interpret the results using outflow models in Section 4, discuss the implications for warm, high-mass superpuffs in Section 5, and conclude in Section 6.

We obtained transit photometry of TOI-1420b with Palomar/WIRC. All events were observed with the beam-shaping diffuser (Stefansson et al. 2017; Vissapragada et al. 2020a), which spreads the light into a top-hat point-spread function with a 3'' FWHM. This allows us to control time-correlated systematics and to integrate for longer on bright sources (improving the duty cycle of the observations). We used the custom metastable helium filter (centered on 1083.3 nm with a bandwidth of 0.64 nm) presented in Vissapragada et al. (2020b).

We observed a full transit of TOI-1420b on UT 2022 June 30 from 05:17 (airmass 1.8) to 12:02 (airmass 1.2) and a partial transit on UT 2022 July 21 from 05:28 (airmass 1.5) to 11:08 (airmass 1.2). The transit duration is T₁₄ = 3.37 ± 0.15 hr (Yoshida et al. 2023). We attempted a partial transit observation on UT 2022 July 28, but poor weather rendered the data unusable for our analysis. All observations were taken with 90 s exposures. To correct the structured background imposed by the narrowband helium filter (further details in Vissapragada et al. 2020b), we obtained 8 dither frames immediately before both transit observations. Additionally, we maintained roughly the same detector positioning for both nights, so we observed the same six reference stars each night (which were selected to be the brightest sources in the field besides TOI-1420, ranging from J = 11.3 to J = 8.8). None of the selected reference stars are known variables, and none exhibited detectable short-term variability in our photometry.

We dark-subtracted, flat-fielded, and background-corrected the data using the methods presented previously in Vissapragada et al. (2020a, 2020b, 2022b) and Greklek-McKeon et al. (2023). We then performed aperture photometry on the target and reference stars using photutils (Bradley et al. 2023). We tested apertures from 5 pixels (125) to 20 pixels (500) in radius. The centers of the apertures were allowed to shift to track variations in telescope pointing. Our guiding was stable on both nights, with centroid variations not exceeding 1.5 pixels. We estimated and corrected any residual local background using annular apertures around each source with inner radii of 25 pixels (625) and outer radii of 50 pixels (1250). The optimal aperture was selected to be the one that minimized the per-point rms scatter in the target star photometry with the average comparison star divided out (removing common-mode systematics) and a moving-median filter applied (removing the transit event). We adopted optimal apertures of 8 pixels (200) for the first night and 7 pixels (175) for the second night.

We first fit each helium curve separately using the exoplanet code (Foreman-Mackey et al. 2021a, 2021b). We fit the target light curve f as the product of a transit light-curve model T and a systematics model S:

$$\begin{eqnarray}&&f=T({R}_{{\rm{p}}}/{R}_{\star },a/{R}_{\star },b,{u}_{1},{u}_{2},P,{t}_{0})\times S.\end{eqnarray} \tag{ 1 }$$

To model the transit, we included a uniform prior on the planet-to-star radius ratio R_p/R_⋆, scaled semimajor axis a/R_⋆, and impact parameter b, uninformative priors (Kipping 2013) on the quadratic limb-darkening coefficients (u₁, u₂), and normal priors on the planetary orbital period P, transit epoch t₀, and stellar radius R_⋆. Normal priors were taken from Yoshida et al. (2023), and are reproduced in Table 1. The planet does not appear to be on an eccentric orbit (e < 0.17 Yoshida et al. 2023) so we fixed the eccentricity to 0. We tracked the depth of the WIRC light curve at mid-transit δ_WIRC numerically rather than using ${({R}_{p}/{R}_{\star })}^{2}$ so that limb darkening was taken into account.

Table 1. Priors and Posteriors for the Light-curve Fits

Parameter	Unit	Prior	Night 1	Night 2	Night 1 + Night 2
P	days	${ \mathcal N }(6.9561063,0.0000017)$	${6.9561063}_{-0.0000013}^{+0.0000013}$	${6.9561065}_{-0.0000013}^{+0.0000013}$	${6.9561064}_{-0.0000013}^{+0.0000013}$
t0	BJDTDB	${ \mathcal N }(2459517.43305,0.00012)$	${2459517.433227}_{-0.000087}^{+0.000087}$	${2459517.433253}_{-0.000087}^{+0.000088}$	${2459517.433238}_{-0.000087}^{+0.000087}$
R⋆	R⊙	${ \mathcal N }(0.923,0.024)$	${0.923}_{-0.024}^{+0.024}$	${0.923}_{-0.023}^{+0.024}$	${0.923}_{-0.023}^{+0.023}$
a/R⋆	⋯	${ \mathcal U }(10,20)$	${16.37}_{-0.36}^{+0.43}$	${16.33}_{-0.36}^{+0.42}$	${16.43}_{-0.38}^{+0.44}$
b	⋯	${ \mathcal U }(0,1)$	${0.429}_{-0.070}^{+0.047}$	${0.429}_{-0.072}^{+0.049}$	${0.414}_{-0.077}^{+0.053}$
Rp /R⋆ (WIRC)	⋯	${ \mathcal U }(0,0.25)$	${0.1393}_{-0.0033}^{+0.0031}$	${0.1477}_{-0.0058}^{+0.0057}$	${0.1409}_{-0.0029}^{+0.0028}$
Rp /R⋆ (TESS)	⋯	${ \mathcal U }(0,0.25)$	${0.1184}_{-0.0017}^{+0.0013}$	${0.1182}_{-0.0018}^{+0.0014}$	${0.1179}_{-0.0018}^{+0.0014}$
u1 (WIRC)	⋯	Kipping (2013)	${0.36}_{-0.23}^{+0.24}$	${0.42}_{-0.24}^{+0.25}$	${0.45}_{-0.22}^{+0.18}$
u2 (WIRC)	⋯	Kipping (2013)	${0.09}_{-0.28}^{+0.34}$	$-{0.05}_{-0.20}^{+0.29}$	$-{0.06}_{-0.19}^{+0.31}$
u1 (TESS)	⋯	Kipping (2013)	${0.44}_{-0.12}^{+0.12}$	${0.41}_{-0.13}^{+0.13}$	${0.40}_{-0.11}^{+0.12}$
u2 (TESS)	⋯	Kipping (2013)	${0.09}_{-0.23}^{+0.25}$	${0.15}_{-0.25}^{+0.27}$	${0.17}_{-0.24}^{+0.25}$
δWIRC	%	⋯	${2.176}_{-0.086}^{+0.088}$	${2.43}_{-0.19}^{+0.19}$	${2.228}_{-0.080}^{+0.081}$
δcomp	%	⋯	${1.568}_{-0.044}^{+0.048}$	${1.553}_{-0.059}^{+0.065}$	${1.556}_{-0.040}^{+0.041}$
δmid	%	⋯	${0.608}_{-0.085}^{+0.085}$	${0.88}_{-0.18}^{+0.18}$	${0.671}_{-0.079}^{+0.077}$

Note. ${ \mathcal N }(a,b)$ indicates a normal distribution with mean a and standard deviation b. ${ \mathcal U }(a,b)$ indicates a uniform distribution between a and b. The WIRC depth δ_WIRC, comparison depth δ_comp, and mid-transit excess depth δ_mid = δ_WIRC − δ_comp are all derived parameters. Posteriors for the detrending weights, WIRC jitter parameters, and TESS error scaling parameters are included in Table 2.

Download table as:  ASCIITypeset image

We included the two-minute cadence TESS light curve of TOI-1420b from Yoshida et al. (2023) in all our fits to provide a reference R_p /R_⋆ for the WIRC light curves. The TESS data also helped us to constrain the degeneracy between b and a/R_⋆ with the precise transit shape. All transit parameters were shared between the WIRC and TESS fits except for R_p /R_⋆ and the quadratic limb-darkening coefficients (u₁, u₂). When we tracked the comparison transit depth δ_comp, we recalculated the TESS light curve with the WIRC limb-darkening coefficients. This excluded the possibility that perceived depth differences were caused by differences in limb darkening between the two bandpasses. The excess depth at mid-transit, which describes the extra absorption due to metastable helium, was then calculated as δ_mid = δ_WIRC − δ_comp.

We modeled the systematics S as a linear combination of detrending vectors s_i including the six comparison starlight curves and the mean-subtracted times (defining a linear trend in time):

$$\begin{eqnarray}&&S=\displaystyle \sum _{i}{w}_{i}{s}_{i}.\end{eqnarray} \tag{ 2 }$$

The weight w_i for each detrending vector s_i in the linear combination was a free parameter with a uniform prior ${ \mathcal U }(-2,2)$. We also tested including the airmass, the centroid offsets, and the water absorption proxy described previously in Paragas et al. (2021) and Vissapragada et al. (2022b) as additional detrending vectors, and calculated the Bayesian information criterion (BIC; Schwarz 1978) to determine whether the increase in log-likelihood justified the addition of these parameters to the model. However, none of these additional detrending vectors improved the BIC, i.e., the marginal increase in final log-likelihood did not justify the addition of any of the parameters, so we left these out of the final models. Finally, we included a photometric jitter term describing the excess white noise in the WIRC data (added in quadrature to the photon noise for each point) and a multiplicative scaling factor for the uncertainties on the TESS data.

In all of the fits, we first optimized toward the maximum a posteriori (MAP) solution with pymc3 (Salvatier et al. 2016). We identified and removed 5σ outliers from the residuals to the MAP fit, iterating the MAP fit and outlier rejection procedure until a stable MAP solution was obtained. Then, we used the No U-Turn Sampler (Hoffman & Gelman 2011) in pymc3 to sample the posterior distributions for the model parameters. Every fit was run with 10,000 tuning steps before taking 10,000 posterior draws (distributed across two chains). To assess convergence, we verified the Gelman–Rubin statistic (Gelman & Rubin 1992) $\hat{R}\ll 1.01$ for all sampled parameters, and we also visually inspected the chains to ensure they were well-mixed. The posteriors of our single-planet fits are given in Table 1.

The transit depths on the two nights were consistent at about the 1σ level. The WIRC transit depth was less precisely constrained on the second night than the first because of the partial phase coverage, but the second night also offered more post-egress baseline than the first night. Given the reasonable agreement in R_p /R_⋆ between both data sets and the complimentary baseline coverage, we decided to fit both WIRC light curves jointly with the TESS light curve. The posteriors of the joint fit are included in Table 1. Figure 2 shows the joint WIRC light curve alongside the phase-folded TESS data. The final WIRC transit was 0.671% ± 0.079% deeper than the TESS transit, indicating substantial extra absorption by metastable helium in the upper atmosphere of TOI-1420b at a confidence of 8.5σ.

Zoom In Zoom Out Reset image size

The rms scatter of the unbinned data was 0.50% on the first night and 0.63% on the second night. These were 1.2 × and 1.3 × the photon noise on each night, respectively, indicating good photometric performance. This is the closest to photon noise we have been able to achieve with Palomar/WIRC using the helium filter (Vissapragada et al. 2022b), likely because six good comparison stars were available in the field. However, when reducing the second night of data, we noticed a large correlated noise component at about a 30 minute timescale in the Allan deviation plots (Figure 3). A correlated feature with this timescale is clearly visible at the end of the WIRC light curve in Figure 2, though the uncertainties in this region of the light curve suggest it is not very significant. This feature was not clearly associated with rapid changes in water absorption proxy, centroid offset, or total flux. Additionally, the second night did not contribute much information to the final excess depth measurement (${\delta }_{\mathrm{mid}}={0.608}_{-0.085}^{+0.085} \%$ for the first night only and ${\delta }_{\mathrm{mid}}={0.671}_{-0.079}^{+0.077} \%$ for the joint fit), so we proceeded without treating this feature in detail.

Zoom In Zoom Out Reset image size

We approached modeling TOI-1420b's helium absorption signal in three stages. First, we estimated the characteristic outflow sound speeds and temperatures that correspond to such a low-density planet. Next, we calculated one-dimensional (1D) models treating the thermodynamics self-consistently rather than assuming an isothermal wind. Finally, informed by the results of the 1D modeling, we calculated a three-dimensional (3D) model that treats the thermodynamics approximately but captures the interaction of the outflow with the tidal and stellar wind environment.

4.1. Characteristic Temperatures

TOI-1420b's low density carries some important implications for the modeling of its outflow. The 1D transonic Parker wind models often adopted for modeling metastable helium absorption (e.g., Oklopčić & Hirata 2018; Lampón et al. 2020; Dos Santos 2023) begin with subsonic flow near the planet that accelerates to supersonic flow at larger radii. A transonic wind requires R_s > R_p , where:

$$\begin{eqnarray}&&{R}_{s}\approx \displaystyle \frac{{{GM}}_{p}}{2{c}_{s}^{2}}\end{eqnarray} \tag{ 3 }$$

(the tidal gravity correction is small for TOI-1420b; see Appendix E of Vissapragada et al. 2022b). The sound speed of the isothermal outflow is:

$$\begin{eqnarray}&&{c}_{s}=\sqrt{\displaystyle \frac{{k}_{{\rm{B}}}{T}_{0}}{\bar{\mu }{m}_{{\rm{H}}}}},\end{eqnarray} \tag{ 4 }$$

where k_B is Boltzmann's constant, T₀ is the isothermal sound speed, $\bar{\mu }$ is the mean molecular weight, and m_H is the mass of a hydrogen atom. The condition R_s > R_p can equivalently be written ${\rm{\Phi }}\gt 2{c}_{s}^{2}$, where Φ = GM_p /R_p is the gravitational potential.

TOI-1420b has a small gravitational potential of Φ ≈ 10¹² erg g⁻¹, so only isothermal outflows with c_s ≲ 8 km s⁻¹ yield transonic velocity structures (i.e., R_s > R_p ). In terms of an isothermal temperature, this corresponds to T₀ ≲ 7700 K for an assumed mean molecular weight of μ = 1 amu. This condition overlaps significantly with the range of isothermal outflow temperatures usually assumed when modeling outflows from hot giant planets (5000 K ≲ T₀ ≲ 15,000 K; e.g., Vissapragada et al. 2022a; Linssen et al. 2022). To avoid the issue of applying the isothermal Parker wind outside the bounds of its validity, we decided to consider 1D models with self-consistent (nonisothermal) thermodynamics and the nature of the transonic ouflows that develop.

4.2. One-dimensional Outflow Models

To self-consistently treat the wind-launching physics, we used the 1D photoionization hydrodynamics code ATES to compute outflow density, velocity, and temperature profiles (Caldiroli et al. 2021). Contrary to an isothermal Parker wind profile, ATES calculates the temperature structure based on the stellar irradiation assuming a pure H/He composition for the wind. With the equilibrium temperature as a boundary condition at R_p , the wind starts hydrostatically and is indeed transonic. We then passed the outflow structure to the NLTE photoionization code Cloudy (Ferland et al. 1998, 2017) to obtain the metastable helium density profile. We finally took Cloudy's metastable helium density profiles and used a custom radiative transfer module (based on Oklopčić & Hirata 2018; Linssen et al. 2022) to calculate the mid-transit spectrum in the stellar rest frame, using the limb-darkening coefficients and transit impact parameter of Table 1. We assumed a hydrogen/helium abundance ratio of 90/10 by number. Finally, we Doppler shift the spectrum from the stellar rest frame into the observatory rest frame, accounting for the system's barycentric radial velocity of −10.33 km s⁻¹ (Yoshida et al. 2023) as well as the Earth's orbital motion relative to the solar system barycenter on the nights of observation. Finally the excess absorption in the WIRC filter passband is calculated following Equation (3) of Vissapragada et al. (2020b).

The XUV flux of the host star is an important ingredient in the simulations, as it governs both the thermal structure in ATES and the metastable helium population in Cloudy. However, the high-energy spectrum of TOI-1420 is unknown (and challenging to measure at d = 200 pc), so we used spectra of similarly active stars as proxies. TOI-1420 is a late G dwarf (T_eff = 5493 K) with $\mathrm{log}(R{{\prime} }_{{HK}})=-4.87$ (Yoshida et al. 2023), indicating similar chromospheric activity to the mean solar value of −4.91 (Mamajek & Hillenbrand 2008). To cover the uncertainty in the stellar XUV spectrum, we repeated the modeling using three different proxies: the MUSCLES spectrum of HD 40307 (an early K star representing "high activity" albeit with a modest $\mathrm{log}(R{{\prime} }_{{HK}})=-4.99$; Mayor et al. 2009; France et al. 2016; Loyd et al. 2016; Youngblood et al. 2016), a quiet solar spectrum (representing "medium activity;" Oklopčić & Hirata 2018; Dos Santos 2023), and the solar spectrum with EUV (<1000 Å) scaled down by a factor of three (representing "low activity").

The high- and medium-activity models predict mass-loss rates of 10^11.25 g s⁻¹ and 10^10.93 g s⁻¹, respectively, but both predict much larger helium signals than we observed (Figure 4). The low-activity spectrum gives the best match to the data at a mass-loss rate of 10^10.58 g s⁻¹, and we therefore adopt this as our preferred XUV model for the star. At the same time, the high- and medium-activity solutions should not be thought of as "ruled out" on the basis of the 1D modeling alone. Some of the absorption in these models comes from high altitudes >4R_p , where shadowing, tidal effects, and interaction with the stellar wind become important for predicting the metastable helium signature (Wang & Dai 2021; MacLeod & Oklopčić 2022). These effects are not treated in ATES, so rather than reducing the activity in the 1D models further, we turned to 3D modeling informed by the results of the self-consistent 1D analysis.

Zoom In Zoom Out Reset image size

4.3. Three-dimensional Outflow Model

In order to assess the impact of the stellar environment on the outflow of TOI-1420b—including the stellar gravity and stellar wind—we perform a 3D hydrodynamic calculation of the escape of a planetary outflow similar to that of the TOI-1420 system. Our model is based on the Athena++ Eulerian hydrodynamic code (Stone et al. 2008, 2020) and the planetary outflow modeling and radiative transfer post-processing software developed by MacLeod & Oklopčić (2022). This simulation does not include magnetic fields, which can suppress the planetary helium signature (Schreyer et al. 2024).

We simulate the interaction of planetary and stellar wind outflows in the corotating reference frame. We assume that the planet and star rotate with the orbital frequency, and adopt the basic system parameters of Table 1. In our model, the planetary and stellar winds are simulated with a nearly isothermal adiabatic index of γ = 1.0001 but have different characteristic temperatures. The planetary outflow is launched from a uniform boundary condition representing the planet and is defined by the hydrodynamic escape parameter ${\lambda }_{{\rm{p}}}\,={{GM}}_{p}/{R}_{p}{c}_{s}^{2}$, where c_s is the sound speed of the isothermal planetary outflow. We adopt λ_p = 5, which corresponds to T ≈ 3200 K for μ = 1 amu. The value of the hydrodynamic escape parameter is chosen to roughly correspond to the 1D ATES model with the 1/3 × solar EUV spectral energy distribution in Figure 4. Similarly, on the basis of the ATES modeling, we adopt the 1/3 × solar EUV as our spectrum for the 3D model spectra.

With these basic parameters of the system set by either the known system properties (e.g., Table 1) or by comparison to the 1D ATES models, we are free to vary the relative strength of the planetary and stellar outflows in our hydrodynamic calculation by specifying the stellar and planetary boundary conditions that generate these flows. This freedom allows us to scale our calculation to reproduce the mid-transit excess absorption. We chose the density of our planetary boundary condition such that the measured emergent outflow rate in steady state is 6.7 × 10¹⁰ g s⁻¹ (10^10.82 g s⁻¹). The stellar boundary condition is set similarly, but with a hydrodynamic escape parameter of λ_* = 10, and a mass-loss rate of 1.6 × 10¹¹ g s⁻¹, approximately 3.2 × the planetary mass-loss rate, or ≈3.5 × 10⁻¹⁵ M_⊙ yr⁻¹. We note that this last choice is somewhat arbitrary, as only much (e.g., order-of-magnitude) stronger or weaker stellar winds are ruled out. Finally, We evolve the system for roughly six orbital periods (41.7 days) and find that the flow reaches a steady state after ∼10 days.

The upper panel of Figure 5 shows the emergent flow in a slice through the orbital plane. The star lies at the origin, and the planet is located in the −x-direction. An inset shows the flow in a region of ±20R_p (roughly ±0.1 au). As the planetary outflow expands into the stellar environment, it is shaped by tidal gravity and stellar wind. Near the planet, outflowing material is initially roughly homogeneous. Outside the planetary Hill sphere of radius ≈4.1R_p , however, the stellar gravity dominates and the nearly radial outflow from the planet is truncated by shocks and sheared into tidal tails leading and trailing the planet. Interaction with the stellar wind also partially shapes these tails and their kinematics in our model, though we note that the value of the stellar wind mass-loss rate is not measured in TOI-1420, so the influence of the stellar wind could be larger or smaller than our nominal assumption (e.g., as described by MacLeod & Oklopčić 2022).

Zoom In Zoom Out Reset image size

We post-process metastable helium level populations and transit spectra in the snapshot shown in Figure 5. For the stellar spectral energy distribution, we adopt the low-activity model described in the preceding section, with 1/3 of the solar EUV flux. We adopt an observer in the general −x-direction but at different viewing angles in order to represent the passage of time as the planet transits. The mid-transit model spectrum is similar to the low-activity model of the preceding section. We then compute the excess absorption in the WIRC filter passband following Equation (3) of Vissapragada et al. (2020b).

The resulting excess absorption light curve is shown in the lower panel of Figure 5. This model achieves similar excess depth at mid-transit to that observed in Figure 2 and shows that the vast majority of absorption is concentrated during the optical transit duration. The extended tails of planetary atmosphere material are effectively hidden by their low optical depth. Thus the region that forms the metastable helium absorption lies within a few planetary radii—near the peak of metastable helium density in the center panel of Figure 4. At these scales, the planetary outflow is nearly spherical and our models thus predict spectra similar to those of the 1D models despite the complexity of the large-scale flow. This optical depth effect can be contrasted with outflows like those observed for HAT-P-32b and HAT-P-67b, where the larger mass-loss rates give the tails greater optical depths, making them readily observable even when compared to the mid-transit absorption (Gully-Santiago et al. 2023; Zhang et al. 2023).

5.1. Can TOI-1420b Retain Its Puffy Hydrogen-rich Atmosphere?

5.2. Is TOI-1420b Tidally Inflated?

While atmospheric retention may not be a problem for high-mass superpuffs, it remains challenging to form planets with such large envelope mass fractions (≳80%; Piaulet et al. 2021; Yoshida et al. 2023) within an in situ core accretion framework (Lee 2019). One possible resolution is that these planets accreted their envelopes far from their stars (≳1 au), where dust-free atmospheres could cool and accrete rapidly (Lee & Chiang 2016; Piaulet et al. 2021). Migration would then be required to bring the planets to their observed positions, and indeed many superpuffs are observed to be parts of resonant chains that suggest such a migration history (Jontof-Hutter 2019). However, TOI-1420b and its doppelgänger WASP-107b do not reside in resonant configurations. WASP-107b is accompanied by a distant companion that might have driven past migration onto the current high-obliquity orbit (Dai & Winn 2017; Piaulet et al. 2021; Rubenzahl et al. 2021), but no additional companions are yet known in the TOI-1420 system (Yoshida et al. 2023).

We therefore turn to tidal inflation as a possible alternative explanation for the low apparent density of TOI-1420b. Although tidal inflation is not required to match the outflow data, the ATES+Cloudy approach used in Section 4 is agnostic to this mechanism. The underlying envelope mass fraction is immaterial to the model so long as the density and temperature at the lower boundary (R_p ) are correctly described. Moreover, the overall escape rate of the flow is relatively insensitive to the boundary density (Salz et al. 2016; Caldiroli et al. 2022). For WASP-107b, the radial velocity data permit a small eccentricity (e < 0.14 at 2σ) and the resulting heating from eccentricity tides could bring the true envelope mass fraction below 50% (Millholland et al. 2020; Piaulet et al. 2021), though this conclusion depends sensitively on the unknown planetary obliquity and tidal quality factor. The radial velocity curve for TOI-1420b permits a similarly small eccentricity (e < 0.17 at 2σ; Yoshida et al. 2023), so the same mechanism may apply.

To investigate the possibility that tidal circularization has heated and inflated the planet, we constructed a model that combines the giant planet thermal evolution models of Thorngren & Fortney (2018) with the equilibrium tides model (e.g., Goldreich & Soter 1966; Mardling & Lin 2002; Jackson et al. 2008). Similarly to Lopez & Fortney (2016), we are not including tidal inspiral or obliquity tides in this model, but we do calculate the amount of heating resulting from circularization to add to the energy balance of the interior. Thus for a constant planetary M_p , bulk metallicity Z_p , and tidal quality factor Q_p , we evolve the eccentricity, semimajor axis, and interior specific entropy over time.

There is a broad space of possible solutions for TOI-1420b, as the model has a number of free parameters (e₀, a₀, Q_p , Z_p , and age) but only strong constraints on a and R_p . Therefore, we provide three illustrative evolution tracks for TOI-1420b in Figure 6. For reference, we show the case in which the planet formed on a circular orbit, which requires an unusually low metallicity of Z_p = 0.13 for a planet of this mass to match the radius (for more details on this nominal model, see Section 4 of Yoshida et al. 2023). This model corresponds to a maximum core mass (assuming all metals are confined to a rocky core) of Z_p M_p = 3.3M_⊕. If instead the planet started with an eccentric orbit, the semimajor axes and eccentricity follow a sigmoid curve with the logarithm of the age. During this process of circularization, a substantial portion of the orbital energy is injected into the planet, resulting in an increase in the planet's radius, which decays as the excess heat is radiated away. The red and yellow curves of Figure 6 show how we may be observing the planet toward the tail end of this process after it has lost most of its initial eccentricity (consistent with the upper limit from Yoshida et al. 2023) but while it still has a moderately hotter interior than it would otherwise. This extra heat allows us to reproduce the observed radius with higher bulk metallicities Z_p = 0.3–0.4, implying core masses that are more in line with core-nucleated accretion (similarly to the case for a tidally heated WASP-107b; Millholland et al. 2020; Piaulet et al. 2021). Different values for the tidal Q_p are correlated with present-day ages in the solution space: having an age over 1 Gyr (which seems likely for this system; Yoshida et al. 2023) would indicate that the planetary Q_p is probably not much less than 10⁶, lest the planet has long finished circularizing.

Zoom In Zoom Out Reset image size

Our calculations demonstrate that tidal inflation is a plausible way to explain the planet's large radius. However, they also demonstrate the strong degeneracies inherent to the problem: all available data for TOI-1420b can be matched equally well in models where the planet is young and circularizing or old and very gas-rich. With relatively uninformative prior information on TOI-1420b's eccentricity (<0.17) and age (<10.7 Gyr Yoshida et al. 2023), these degeneracies cannot be broken. Future efforts to more precisely constrain the eccentricity and/or age of the planet would therefore be highly valuable for understanding TOI-1420b's evolutionary state.

5.3. Can High-altitude Hazes "Deflate" TOI-1420b?

The flat transmission spectra of lower-mass superpuffs are well-modeled by high-altitude aerosols, which may also be responsible for their large apparent radii (Kawashima et al. 2019; Gao & Zhang 2020; Ohno & Tanaka 2021). It is worth considering this hypothesis in more detail for the higher-mass superpuff TOI-1420b, as the model requires strong outflows (similar to the one we have detected) to carry aerosols to high altitudes. While high-altitude aerosols are ruled out as a possibility for WASP-107b (molecular features are detected lower in its atmosphere; Kreidberg et al. 2018; Spake et al. 2018; Dyrek et al. 2024), transmission spectroscopy sensitive to the lower atmosphere of TOI-1420b is not yet available.

Hazes can persist at μbar to nbar pressures in atmospheres undergoing outflows (Wang & Dai 2019; Gao & Zhang 2020; Ohno & Tanaka 2021), and so the 1 nbar radius should represent the maximum extent to which hazes can increase a planet's radius from a clear atmosphere case. As such, to test whether TOI-1420b could be a smaller planet with optically thick high-altitude hazes, we compute the 1 nbar radius of this planet using the atmospheric and planetary structure model of Gao & Zhang (2020). Briefly, the model assumes a core with an Earth-like composition and a mass–radius relation adopted from Zeng et al. (2019), which lies beneath a gas envelope of solar composition and a simplified thermal structure consisting of a deep convective zone and an upper radiative zone that is connected to the convective zone at the radiative–convective boundary (RCB). The convective zone structure is solved following Owen & Wu (2017) assuming hydrostatic equilibrium and a polytrope equation of state with a power of 7/5 to represent a largely H₂ composition. At the RCB, the temperature–pressure profile transitions smoothly to a radiative profile following Piso & Youdin (2014), with the temperature of the RCB set to the equilibrium temperature of TOI-1420b. For the opacity in the radiative zone, we use the tables from Freedman et al. (2014). We refer the reader to Gao & Zhang (2020) for more details regarding the model.

In Figure 7, we show that the 1 nbar radius of TOI-1420b does not reach the observed radius for any atmospheric mass fractions <30% and intrinsic temperatures <100 K, even when the uncertainty in the mass is taken into account. These limits roughly correspond to those of Neptune-mass planets with the age and insolation of TOI-1420b (Helled et al. 2011; Lopez & Fortney 2014). Therefore, it is unlikely that TOI-1420b is a smaller planet that has been made to look larger with a high-altitude haze layer.

Zoom In Zoom Out Reset image size

Many lower-mass superpuffs (M_p < 10M_⊕) are expected to have atmospheric lifetimes shorter than their ages, while higher-mass (M_p > 10M_⊕) superpuffs are expected to be stable to catastrophic atmospheric loss even when strongly irradiated. To test this expectation, we searched for an outflow in the atmosphere of TOI-1420b, a warm (T_eq = 960 K), higher-mass (M_p = 25.1M_⊕) superpuff with a density of just ρ = 0.08 g cm⁻³ (Yoshida et al. 2023). We observed one full transit and one partial transit with the narrowband helium filter on Palomar/WIRC and compared the resulting light curve to the TESS light curve. We found that the transit in the helium bandpass was consistently 0.671% ± 0.079% deeper than the TESS transit, indicating excess absorption by outflowing metastable helium at a confidence of 8.5σ.

We first interpreted our observations using the 1D photoionization hydrodynamics code ATES coupled with Cloudy, which allowed for a self-consistent treatment of the outflow thermodynamics. The model that best matched our data used a low-activity proxy for the host star's XUV spectrum (1/3 × solar XUV) with a hydrodynamic escape parameter λ_p ≈ 5. Using this XUV spectrum and escape parameter, we then calculated a 3D isothermal outflow model using Athena++. This model gave a good match to the observed excess absorption at mid-transit with a mass-loss rate of $\dot{M}=4.9\,\times {10}^{10}$ g s⁻¹. The corresponding mass-loss timescale is ${M}_{p}/\dot{M}\,\approx 70\,\mathrm{Gyr};$ despite its remarkably low density, the mass-loss rate of TOI-1420b is quite modest. The situation is similar for the other warm high-mass superpuff WASP-107b, which has ${M}_{p}/\dot{M}\approx 30\,\mathrm{Gyr}$ (Wang & Dai 2021). Our observations confirm that superpuffs with M_p > 10M_⊕ are easily able to retain their massive atmospheres over long timescales. These planets may therefore have formed with unusually massive envelopes, likely accreting gas outside the water ice line before migrating to their present positions (Stevenson 1984; Venturini et al. 2015; Lee & Chiang 2016; Piaulet et al. 2021).

We also explored a number of alternative hypotheses for TOI-1420b's large radius. We presented some illustrative evolution calculations showing that all available data for TOI-1420b could be matched at a larger bulk metallicity (more consistent with models of core accretion in the inner disk; Lee 2019) provided the planet recently underwent tidal circularization. On the other hand, we found that high-altitude hazes are disfavored as an explanation for TOI-1420b's large radius. The planet is simply too massive, resulting in a relatively small scale height among the superpuffs; even optically thick hazes at nbar pressure for a Neptune-like planet would fall well short of explaining TOI-1420b's large radius. This bodes well for future atmospheric characterization efforts: we do not anticipate that atmospheric features will be completely obscured by hazes, similar to the case for WASP-107b (Kreidberg et al. 2018; Spake et al. 2018; Dyrek et al. 2024).

TOI-1420b's helium absorption is the largest signal we have observed at high confidence with Palomar/WIRC thus far (Vissapragada et al. 2022b). We encourage other investigators to follow up on the helium absorption for this planet with high-resolution spectroscopy, which would provide additional constraints on atmospheric modeling. With such a strong signal this should also be a good candidate to search for a Doppler-shifted line core and/or extended egress, which is indicative of an outflow morphology sculpted by stellar winds in WASP-69b and WASP-107b (Nortmann et al. 2018; Spake et al. 2021; Wang & Dai 2021; MacLeod & Oklopčić 2022; Tyler et al. 2024). Our light curve of TOI-1420b appears symmetric and is thus well-fit by a 3D model where the stellar wind is not strong enough to shepherd the outflowing material into a comet-like tail, but tail features can be challenging to capture at low effective resolving power from the ground (Vissapragada et al. 2020b; Paragas et al. 2021; Fu et al. 2022). High-resolution studies of TOI-1420b would therefore enable more direct morphological comparisons to planets like WASP-107b, and may even enable comparative studies of the stellar wind strengths in these systems.

Lower-mass superpuffs are expected to be even more susceptible to atmospheric escape than TOI-1420b and WASP-107b. In particular, there are a number of planets in Figure 1 that are similar in equilibrium temperature to WASP-107b, including TOI-216b, Kepler-79d, Kepler-18d, and Kepler-223e. If their host stars are active, these planets might also exhibit triplet helium absorption. While many superpuffs reside in systems that are too faint for current facilities, a number of suitable targets have been discovered in bright systems, including K2-24 c, HIP 41378f, and TOI-216b (Petigura et al. 2018; Santerne et al. 2019; McKee & Montet 2023). Although the transit durations for these planets tend to be long, helium absorption could still be constrained with multinight ground-based campaigns (Gully-Santiago et al. 2023; Zhang et al. 2023) or space-based observations with JWST (Fu et al. 2022; Dos Santos et al. 2023). By characterizing the upper atmospheres of lower-mass superpuffs, we can gain greater insight into the nature of these remarkably low-density worlds.

We thank Paul Nied and Diana Roderick for their assistance with telescope operations. We additionally thank Fei Dai, Ryan Rubenzahl, and Jeremy Drake for helpful conversations. This work was supported by telescope time allocated to NASA-NSF Exoplanet Observational Research (NN-EXPLORE) partnership through the scientific partnership of the National Aeronautics and Space Administration, the National Science Foundation, and the NOIRLab. This research made use of photutils, an astropy package for detection and photometry of astronomical sources (Bradley et al. 2023). This research made use of exoplanet (Foreman-Mackey et al. 2021a, 2021b) and its dependencies (Astropy Collaboration et al. 2013; Kipping 2013; Salvatier et al. 2016; Theano Development Team 2016; Astropy Collaboration et al. 2018; Kumar et al. 2019; Luger et al. 2019; Agol et al. 2020).

Facility: ADS - , NASA Exoplanet Archive - , Hale 200 inch (WIRC) - .

Software: photutils (Bradley et al. 2023), exoplanet (Foreman-Mackey et al. 2021a, 2021b), pymc3 (Salvatier et al. 2016), arviz (Kumar et al. 2019), numpy (Harris et al. 2020), scipy (Virtanen et al. 2020), astropy (Astropy Collaboration et al. 2013, 2018, 2022), p-winds (Dos Santos 2023), ATES (Caldiroli et al. 2021), Cloudy (Ferland et al. 1998, 2017), Athena++ (Stone et al. 2008, 2020).

Here, we provide a corner plot for the key fitting parameters of the final joint fit in Figure 8 and a table summarizing the posteriors for detrending weights and error scaling parameters in Table 2.

Zoom In Zoom Out Reset image size

Table 2. Priors and Posteriors for the TESS error Scaling Parameter, Jitters, and Detrending Weights

Parameter	Prior	Night 1	Night 2	Night 1 + Night 2
kTESS	${ \mathcal U }(0.5,1.5)$	${0.9158}_{-0.0022}^{+0.0021}$	${0.9158}_{-0.0022}^{+0.0022}$	${0.9158}_{-0.0022}^{+0.0021}$
σWIRC, 1	${ \mathcal U }({10}^{-6},{10}^{-2})$	${0.00204}_{-0.00067}^{+0.00053}$	⋯	${0.00202}_{-0.00065}^{+0.00053}$
w1,1	${ \mathcal U }(-2,2)$	${0.073}_{-0.035}^{+0.034}$	⋯	${0.070}_{-0.034}^{+0.035}$
w2,1	${ \mathcal U }(-2,2)$	${0.719}_{-0.078}^{+0.079}$	⋯	${0.719}_{-0.075}^{+0.076}$
w3,1	${ \mathcal U }(-2,2)$	${0.0072}_{-0.0393}^{+0.0402}$	⋯	${0.0083}_{-0.0414}^{+0.0410}$
w4,1	${ \mathcal U }(-2,2)$	${0.014}_{-0.036}^{+0.037}$	⋯	${0.012}_{-0.036}^{+0.037}$
w5,1	${ \mathcal U }(-2,2)$	${0.129}_{-0.051}^{+0.053}$	⋯	${0.138}_{-0.051}^{+0.049}$
w6,1	${ \mathcal U }(-2,2)$	${0.072}_{-0.038}^{+0.038}$	⋯	${0.068}_{-0.039}^{+0.039}$
w7,1	${ \mathcal U }(-2,2)$	${0.0146}_{-0.0052}^{+0.0051}$	⋯	${0.0152}_{-0.0052}^{+0.0050}$
σWIRC, 2	${ \mathcal U }({10}^{-6},{10}^{-2})$	⋯	${0.00361}_{-0.00049}^{+0.00051}$	${0.00361}_{-0.00049}^{+0.00050}$
w1,2	${ \mathcal U }(-2,2)$	⋯	${0.254}_{-0.092}^{+0.095}$	${0.250}_{-0.094}^{+0.093}$
w2,2	${ \mathcal U }(-2,2)$	⋯	${0.44}_{-0.12}^{+0.13}$	${0.43}_{-0.12}^{+0.12}$
w3,2	${ \mathcal U }(-2,2)$	⋯	${0.093}_{-0.047}^{+0.046}$	${0.098}_{-0.045}^{+0.044}$
w4,2	${ \mathcal U }(-2,2)$	⋯	${0.083}_{-0.050}^{+0.050}$	${0.093}_{-0.051}^{+0.050}$
w5,2	${ \mathcal U }(-2,2)$	⋯	${0.034}_{-0.090}^{+0.087}$	${0.027}_{-0.087}^{+0.090}$
w6,2	${ \mathcal U }(-2,2)$	⋯	${0.121}_{-0.060}^{+0.060}$	${0.116}_{-0.061}^{+0.062}$
w7,2	${ \mathcal U }(-2,2)$	⋯	$-{0.006}_{-0.016}^{+0.016}$	${0.007}_{-0.013}^{+0.013}$

Download table as:  ASCIITypeset image